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The strong-coupling perturbation theory ol the Hubbard model is presented and carried out 
' to order (t/U)^ for the one-particle Green function in arbitrary dimension. The spectral weight 

j4(k,tj) is expressed as a Jacobi continued fraction and compared with new Monte-Carlo data 
of the one-dimensional, half-filled Hubbard model. Different regimes (insulator, conductor and 
short-range antiferromagnet) are identified in the temperature-hopping integral (T, t) plane. This 
work completes a first paper on the subject (Phys. Rev. Lett. 80, 5389 (1998)) by providing 
details on diagrammatic rules and higher-order results. In addition, the non half-filled case, infinite 
resummations of diagrams and the double occupancy are discussed. Various tests of the method are 
also presented. 
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I. INTRODUCTION 



The study of strongly-correlated electrons has become in the last decade one of the most active fields of condensed 
matter physics. The electronic properties of an increasing body of maierials cannot be described adequately by Lan- 
dau's theory of weakly- interacting quasiparticles (Fermi liquid theory) .El'B Best known are the high-Tc superconductors 
Q ' and organic conductors. In both cases, a strong anisotropy and a narrow conduction band contribute to make the 
O i effects of interactions between electrons (mainly Coulomb repulsion) dramatic. 

From a theoretical point of view, quasi-one dimensional (e.g. organic conductors) and quasi-two dimensional (e.g. 
_^ high-Tc superconductors) systems are quite different. In one dimcasion, it is possible to solve satisfactori|j|r a great 



^ , number of models : lattice Hamiltonians such as the Hubbard models 
fS| ' nonperturbative results can be obtaiaed for the Tomonaga-Luttingerl 



fm be solved exactly by Bethe Ansatz, 3 whereas 
and the (7-ology models from bosonizationQ and 
renornpJigation-group techniques B'EI Conformal field theory has been applied as well, in particular to the Hubbard 
model,ll3ll2l whose Bethe Ansatz solution has limited practical utility. A unified phenomenology of the -sp-called 
^"^^ Luttinger liquids emerges from these works, whose most striking feature is probably spin-charge separationJlj 

On the other hand, the above methods and results cannot be generalized to the case of two dimensions, relevant to 
the Cu02 planes present in all high-Tc superconductors. With a half-filled conduction band, all the parent compounds 
, of the cuprates are antiferromagnetic insulators, and the main challenge is to understand, first towards which kind of 
d ' metal they evolve upon doping, and secondly whether superconductivity can occur, without a phonon-mediated coifc 
^ ' pling, in the vicinity of this antiferromagnetic phase. Moreover, the linear temperature dependence of the resistivityt^ 
I ' is a strong experimental indication that the normal phase of high-Tc superconductors is not a Fermi liquid. 
' O ' In this context, the Hubbard modeH 

a : 

3 : n^-tY^clc.^+UY, 4t 4l C^l C^^ (1) 



has spurred renewed interest for its ability to account for antiferromagnetic correlations and the Mott metal-insulator 
transition. It is also the simplest model of interacting electrons. In Eq. (|^), t represents the hopping amplitude 
between two neighboring sites in a tight-binding approximation, and U the strength of the very effectively screened 
(and thus taken to be local) Coulomb repulsion. At half-filling and for t U, the Hubbard model is equivalent to 
a Heisenberg model with antiferromagnetic exchange J = At^ /U , and its ground state has long-range Neel order in 
any dimension d > 2. A Mott transition can occur towards a metallic state either upon doping, or by increasing the 
ratio t/U . Thus, the Hubbard model is the prototype of strongly-correlated electrons systems and has beep|intensely 
studied, although a more realistic model of the Cu02 planes of the cuprates should involve several bands.li3 

Solving the Hubbard model is a difficult problem by itself, even in dimension = 1, where the complexity of 
the Bethe Ansatz solution prevents one from actually computing piost physical quantities. For example, its suectrjal 
weight is knmjiruanly in the limit U 00 and at zero temperature. Ej In any dimension d > 2, only apprcajimateB^E^I E£l 
LlE2rc3 methods are available. Only in the limit d — > 00 do important simplifications occur,E3 allowing to 
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compute most physically interesting quantities in an essentiallv exact wayB3'^ In particular, the Mott transition has 
been studied in great detail, and is still under discussion E3'L2I Strong coupling perturbation theory, which considers 
the interaction term of the Hamiltonian as dominant, and the kinetic term as a perturbation, has beea_spmewhat 
neglected so far. Perturbative series for the thermodynamical potential have been obtained up to order t^'jETObut do 
not yield dynamical correlation functions. An original method effectively achieving infinite summations of terms, tha 
so-called Grassaaannian Hubbard-Stratonovich transformation, was introduced independently by C. BourbonnaisEfl 
and S. Sarker,L3 but contained great difhculties which were overcome only recentlyE3 by the authors. The purpose 
of the present paper is to develop the ideas of Ref. ^ in greater detail, and to present new results that have been 
obtained in the meantime. 

Although we are naturally interested in the two-particle correlation functions (magnetic susceptibility, conductivity, 
compressibility) and the phase transitions of the Hubbard model, it has proven simpler to first leijtiGjdate one-particle 
properties. Furthermore, spin-charge separation in d = 1 is mpatly visible in the spectral weight,E3E3 and it is possible 
to gain significant insight into the metal-insulator transitions and antiferromagnetic correlations from the spectral 
weight alone. Thus, we will focus on the one-particle properties of the Hubbard model throughout most of this paper. 
From an experimental point of view, angle-resolved photoemission spectroscopy (ARPES) is the probe of chpijce to 
measure_the spectral weight. Experiments have receatly been conducted on quasi one-dimensional SrCu02E3 and 
NaV205L2l, and on quasi two-diaiensional SrCuCl202cil antiferromagnetic compounds, and further analyzed with the 
help of exact diagonalizations.L3 

In Sect. |l|, we present the strong-coupling expansion for dynamical correlation functions, and provide explicit 
examples of diagram calculation in Appendix]^. The method is applied to the half- filled Hubbard model in Sect. HI, 
and the resulting spectral weight and double-occupation are compared to Monte-Carlo data in Sect. IV. Partially 
self-consistent solutions, involving an infinite sum of diagrams, are investigated in Sect. and the question of doping 
addressed in section VI. Appendix]^ provides the atomic one- and two-particle correlation functions, necessary to 
apply the method, and appendix O presents a practical test of the method on a toy model. Higher-order terms of the 
expansion are given in appendix O. 



II. THE STRONG-COUPLING EXPANSION 



In this section we derive the strong-coupling expansion of correlation functions for a wide class of Hamiltonians. 
We first specify the form that the Hamiltonian should have in order for the method to work. Then we introduce the 
Grassmannian Hubbard-Stratonovich transformation, and the diagrammatic perturbation theory itself. 



A. The Hamiltonian 



Consider a Hamiltonian Ti. = TiP + H^, where the unperturbed part is diagonal in a certain variable i (for 
instance a site variable), and let us denote collectively by a all the other variables of the problem (for instance a spin 
variable). From now on we will call i the "site variable" and a the "spin variable" for definiteness, though they may 
represent any set of quantum numbers. This Hamiltonian describes the behavior of fermions, and we suppose that 
it is normal-ordered in terms of the annihilation and creation operators c^jj . Ti!^ may be written as a sum over i of 
on-site Hamiltonians involving only the operators c^P at site i: 

H"-5]/i.(4,c..). (2) 

i 

Whenever doing actual calculations, we will use the Hubbard model. For the latter, Ti^ corresponds to the atomic 
limit, that is 

ht{c\^,c^a) = C/cj^4^Q|Q| . (3) 

We will use the notation u = U/2 throughout this paper for convenience. We suppose that the perturbation Ti^ is a 
one-body operator: 

w'-EE^^^^c,., (4) 

a ij 
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with V a Hermitian matrix. Here, we suppose in addition that the perturbation is diagonal in the spin variable, but 
this needs not be the case. For the Hubbard model, Ti.^ represents the kinetic term, and V is the matrix of hopping 
amplitudes. 

Introducing the imaginary-time dependent Grassmann fields 7iaiT),"f*^{T), the partition function at some temper- 
ature T = may be written using the Feynman path-integral formalism: 



(5) 



In order to lighten the notation, we will use the first few latin letters (a, 6, ...) to denote sets such as (i, ct, r), and use 
bra-ket notation. For instance: 

r dT J2 Vrnl{r)lja (r) = ^-"^alb = {iMl) ■ (6) 

ija ab 



B. The Grassmannian Hubbard- Stratonovich transformation 



At this point, one could use standard perturbation theory and expand the ^-matrix in terms of unperturbed 
correlation functions. Due to the absence of Wick theorem in the case of a nonquadratic unperturbed Hamiltonian, 
this approach does not lend itself to a satisfactory diagrammatic theory. One cannot define one-particle irreducibility 
for the diagrams_and one has to be very carefuJ in order to avoid double counting of certain contributions. For example. 
Pan and Wangj23 and Bartkowiak and ChaoEj had to deal with over five hundred different diagrams in order to obtain, 
the thermodynamic potential of the Hubbard model up to fourth order. These problems were solved by Metzner,a 
who organized the perturbation series as a cumulant expansion, and formulated diagrammatic rules with unrestricted 
sums over momenta. In this subsection we show that Metzner's results can be obtained in a straightforward fashion 
and even further simplified byi-^eans of a simple transformation on the partition function. This-|transformation was 
first proposed by BourbonnaisEJ and applied (at zeroth order) to the Hubbard model by Saiifir.L^ D. Boies et al. also 
used it to study the stability of several Luttinger liquids coupled by an interchain hopping.C3 

The Grassmannian Hubbard-Stratonovich transformation amounts to expressing the perturbation part of the ex- 
ponential as the result of a Gaussian integral over auxiliary Grassmann fields ^ia-{T),^*^{T) as follows: 



[d^^'diP] e('^l^ '\i') + W-f)+iiW = det(y^i) e-<^l^l'^> , (7) 
In terms of the auxiliary field, the partition function becomes, up to a normalization factor: 

Z = ZoJ [dV*rfV']e<'^l^"'l'^)(eW'l'^>+<'^l'^'>^^ , (8) 

where (...)o means averaging with respect to the unperturbed Hamiltonian. Denoting by (...)o ^ the cumulant averages, 
and owing to the block-diagonality of 7i°, the average in Eq. (^) can be rewritten as: 

CO ^ 13 R. 

'^^pEt^ E / ]l'^'^ld'^l'^Ui'^l)■■^ta,,irR)i^^a'JT'Ji)..iJ,„,^{T^ (9) 
R=l ^ i{ai,a[}''^ 1 = 1 

We will denote by G^f'^ln = (~)^(7ai ■■7afi76 ■■Ih )n i ^ ^^'^ various Green functions of the unperturbed Hamiltonian. 



/0,(c) 

The partition function now takes the familiar form 



Z cx y [d^di^] exp J -Soir, V'] - E ^int[^^ ^] I ' 



(10) 



R=l 

where the action has a free (Gaussian) part 

So[r,^]^~{i'\v-'\^p) , (11) 



3 



and an infinite number of interaction term; 



^int[V'*,^] = E Va,..^a,K-<Gfi-H ■ (12) 

The primed summation in Eq. (|l2|) reminds us that all the fields in a given term of the sum share the same value of 
the site index. The free propagator for the auxiliary fermions is just V, and we may now use Wick's theorem to derive 
diagrammatic rules in order to treat the interaction terms perturbatively. Appendix |^ gives a thorough description 
of the diagrammatic rules as well as two explicit examples of application to the one-particle Green function and the 
thermodynamical potential. 

One may wonder why we did not include the first interaction term into the free part of the action 5*0 [V'*, since 
it is quadratic in the field tp. The reason is that we want to count precisely the order of a given diagram in the 
perturbation V. By separating the free and interacting parts of the action as in Eqs. ( pT| , p^ , the order of a given 
diagram is simply the number of ip propagators, whereas the alternate method mixes all powers of V from zero to 
infinity. This question will be discussed more thoroughly in Sect. M below. 



C. Electron Green functions 



We now show how to deduce the correlation functions of the original electrons (7) from those of the auxiliary fermions 
(rp) to which the perturbation theory applies. When coupling the electron field 7 to an external Grassmannian source 
J*, Ja, the partition function takes the form (very similar to Eq. (jsl)): 



z(J^J) 

A 2_R-point correlation function reads 



Noticing that 



(-)^(7ai-7a„7b„-7?i) 



"(S(e<''^+'^I'>'>+<^I'''+"'>)Q 



SJ*,..SJ*JJb^..6J, 



,(V|v-i|v) 



:.{i>+Jh)+{'(\->P+J) 



.SJ, 



bi 



6{, 



.7 = 
.J*=0 



we perform 2R integrations by parts, and obtain: 



/o 



(13) 



(14) 



(15) 



(16) 



which expresses in terms of correlation functions of the tp field. For instance, a straightforward calculation gives 
the relation between the Green functions Qab = ~(7a7b) snd Vab — "(V'aV'b): which we write in matrix form 

g ^ ~v-^ + v-^vv-^ . (17) 

If r denotes the self-energy of the auxiliary field, then 

g^{T-^-vy\ (18) 

Likewise, a connected 2-point correlation function of the original fermions is the corresponding amputated, connected 
2-point correlation function of the auxiliary fermions. 

qWc _ ^^r~'i-\ -iillc 



'aia2,6ift2 



(19) 



where a summation over repeated indices is implicit on the right-hand side. 

To conclude this section, let us summarize our line of thought up to now. We wish to compute physical quantities 
within a strongly correlated fermion model (say, the Hubbard model). Given that the energy of the interaction is 
greater than the bandwidth, we want to build a strong-coupling expansion. Ordinary perturbation theory turns 
out to be quite cumbersome, but the simple transformation (^ restores Wick's theorem and a (nearly standard) 
diagrammatic approach for the auxiliary fermions. Finally, the simple relations ( pT|jT^Jl9| ) make the connection back 
to electron Green functions. 
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III. APPLICATION TO THE HALF-FILLED HUBBARD MODEL 



We apply the method presented in the previous section to the Hubbard model at half filling in any dimension. 
When doing so, we first have to deal with an unexpected causality problem, whose solution can be viewed as part of 
the method itself. Then we compute the spectral function and the density of states, and discuss their relevance to 
the Mott metal-insulator transition and to antiferromagnetic correlations. 



A. The one-particle Green function up to order three 



The simplest dynamical quantity amenable to our method is the one-particle spectral function: 

A(k,uj)= lim -2lmg(k,uj + if]), (20) 

It is the simplest energy and wavevector resolved correlation function. As such, it contains a lot of physical information, 
far more detailed than the density of states, which is only its momentum-integrated version. It is also the best tool 
to investigate the Fermi liquid character of the system, i.e., whether or not it is dominated by a narrow peak at the 
Fermi level. Finally, it is very badly known from a theoretical (and numerical) point of view, whereas experiments of 
angle-resolved photoemission have recently become more reliable and precise. To obtain the Green function, we have 
to compute the self-energy T of the auxiliary fermions. From now on we will use the nearest-neighbor Hubbard model 

'H = -tY^ cl^Cja + 2u^4^cT^Qxc,| (u=^\ (21) 



which means that V^(k) = — 2te(k), where c(k) = J2m=i coskm- Throughout this section, we work at half-filling by 
setting the chemical potential to fj, = u = U/2. The simplest diagram contributing to F, of order zero in t, is just the 
atomic Green function 

F(0) (k, luj) = Gituj) = .. , , (22) 

where ilo stand for a complex frequency (we keep the symbol lo for a real frequency). This leads to the following 
approximate Green function: 

Q^'\^^) = -—2 \ • (23) 



{iuj) - u 



2te(k) 



One recognizes in Eq. ( p3D the result of the Hubbard-I approximation, whose properties are well known. The original 
atomic levels at ±u are spread out by the hop ping term into two symmetric bands called Hubbard bands. The lower 
band has band edges at ±2td — ^ {2td)^ + v? and never reaches the Fermi level, however large t may be. Thus, this 
approximation fails to describe the Mott metal-insulator transition. But in the method used here, the Hubbard-I 
approximation simply stems from the zeroth order value of F, and we know a systematic way of improving it. 



o 




FIG. 1. Diagrams contributing to the self-energy F of the auxiliary fermions up to order . 

The diagrams contributing to F up to order are presented in Fig. |l]. Actually, there are more diagrams at this 
order, but they vanish due to the precise form of T^(k). To compute these diagrams, we need the explicit expression of 
^^(71(72 ct30-4(*'^1j *'^2; ^1^3) jw4). The calculation of this atomic correlation function is a bit cumbersome, though simple 
in principle. The main result as well as several useful remarks are presented in Appendix |b| For the time being, let 
us focus on the result: 
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By injecting the above approximate value of F into Eq. ( p^ for the Green function, one obtains a rational function of 
iw, presenting a finite number of poles in the complex plane. Such a result is already somewhat disappointing since 
it cannot account for a continuous spectral weight. But more seriously, the Green function presents pairs of complex 
conjugate poles away from the real axis, and as a consequence does not verify Kramers-Kronig relations. Furthermore, 
for a wavevector verifying F(k) = 0, the Green function is equal to F itself, and has high-order poles at ±u, leading 
to negative spectral weight. So by simply replacing F by its approximate value to third order in t in Eq. (^8|), one 
ends up with a noncausal Green function, whose spectral function is neither positive, nor normalized. 

Any Green function acceptable on physical grounds should be causal and have a positive spectral weight, i.e., be 
a sum (finite or infinite) of simple real poles with positive residues. We call such a function Lehmann representable 
since these properties are obvious when considering the Lehmann representation of the Green function. The way out 
of the problem is to seek a Lehmann representable function resembling the approximate expression of Q as much as 
possible. Since we know Q up to third order only, any function that differs from it by terms of order is a priori as 
good an approximation. A systematic way of building such a function is provided by a theorem reported in Ref. 
a rational function is Lehmann representable if and only if it can be written as a Jacobi continued fraction 

Gj{i(^) = — -T — — -T — ■■• ——J- , (25) 

lUJ + bi— lUJ + 02— lUJ + On 

with 



hi real and a; > 0. (26) 

In any finite system, the exact Green function can be cast into the form of Eq. (^5|) with a finite value of n (the number 
of floors of the continued fraction). If we consider the coefficients a; and hi as functions of t and expand Qj{iuj) in 
powers of t, we destroy its continued fraction structure, and may well obtain an unacceptable approximation. On the 
other hand, if we expand the coefficients themselves to some finite order in t and leave them where they belong in 
Qj{iuj), we can still expect their truncated Taylor series to verify conditions Eq. ([2^), and consequently the approximate 
Gj{iuj) obtained this way to be Lehmann representable. In addition, their exists a double recurrence relation giving 
the ai and hi starting from the moments of the function. Since we know the exact Taylor expansion of G up to order 
t^, we know all its spectral moments to the same precision. So all we have to do is to work out the recurrence relations 
and compute the coefficients of the continued fraction up to the best precision available. Once an ai is found to be 
zero to the best of our knowledge - that is, if it contributes to the full Gj, given the preceding coefficients, at an order 
higher than the working precision {t^ in this section) - we have to truncate the continued fraction. In all the cases 
that we treated, such an ai always occurred rather quickly (the deepest continued fraction we had to deal with had 
eight fioors). This means that the density of poles is always weak, and the poles remain essentially well separated: 
they cannot account for the precise shape of a continuous spectral weight, but can sample it adequately (the moments 
of the distribution are well represented). More important, the approximate Green function thus obtained is causal. 

In order to test both the diagrammatic theory and the continued fraction representation, we treated the exactly 
solvable case of the atomic limit itself, away from half-filling. Indeed, setting the chemical potential to /i = m + tg, 
one can either compute the exact Green function, then expand the solution in powers of to i or notice that the shift of 
chemical potential has the form of a zero-range hopping and obtain the to expansion from diagrams. This procedure is 
presented in Appendix and leads to the following main conclusions: the diagrams actually give the correct answer, 
which however presents the same causality problems as already mentioned, and the continued fraction representation 
properly builds back the poles and the residues of the exact solution. We also tested our approach on the (again 
exactly solvable) two-site problem, where the same conclusions apply. 

The above considerations also shed new light on the usual weak-coupling perturbation theory. In that case too, one 
gets multiple poles when truncating the series for G, and the way out of this difficulty is to use Dyson's equation - valid 
because Wick's theorem applies directly to the original fermions - and to compute the self-energy. If the self-energy 
is Lehmann representable, that is, if it possesses an underlying continued fraction structure, then the Green function 
will inherit this structure from it and will be Lehmann representable too. We did not find any definite answer as to 
why the self-energy itself, as obtained from a few diagrams, turns out to be acceptable in general, but there are some 
plausible explanations. One is that the unperturbed case already presents a continuum of levels, instead of only two. 
So if a double pole appears in the expression of a diagram, it is likely that its isolated contribution will be negligible in 
the thermodynamic limit. Also, the vertices of weak-coupling theories are often local in time or, if retarded, depend 
only on two times, whereas here the vertices have a full dependence on all frequencies entering them, which is quite 
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peculiar. Finally, nothing proves that a finite-order weak-coupling self-energy is always Lehmann representable, and 
counter-examples may well exist. 

Applying the above procedure to the result ( p4[ ) yields the following continued fraction: 

_ 1 6t^d V? 

•'^ ' ~ iw-K 2te(k)- - 3/3*3 tanh (/3M/2)c(k) /u- iu ^ 2tc{W, j d- iuj + tcCk)/d ' ^ 



which has exactly the same Taylor expansion as the exact Green function up to order t^, verifies the conditions ([2^), 
and is normalized to unity. The next two subsections show that the spectral weight deduced from Eq. ( p^ describes 
the Mott metal-insulator transition, and the appearance of strong antiferromagnetic correlations at low temperature 
in the Hubbard model. 



B. The Mott transition 



Strictly speaking, there is no rigorous definition of the Mott transition in terms of one-particle properties only. But 
one can use as a heuristic criterion (even at nonzero temperature) the appearance of spectral weight at the Fermi 
level. With the normalization (20) of the spectral weight, the density of states has the following expression: 



N{uj) 



d'^k 



A(k,c.). 



(28) 



We observe that in the density of states, as t increases from zero, the two symmetric Hubbard bands - reduced to 
two delta functions at ±it in the atomic limit - widen, and eventually mix for t beyond some critical value. At high 
temperature {T > u), the gap is closed by excitations of momentum k = (0, 0) and k = (vr, tt), making it possible 
to compute the exact value of the critical hopping: 



(29) 



This gives Uc — 3.2i for c? = 1, to be compared with Uc — 3.5t found in the Hubbard-III approximation. Note that tc 
scales properly with the dimension, and leads to-LL_= l.SGt* (with t* = 2t^/d) in the limit of infinite dimension. This 
critical value of interaction strength is too largeuZrED for two reasons. First, when rf — > cx) our criterion corresponds to 
subbands that meet with an exponentially small density of states, therefore not yet truly closing the gap. Secondly, 
when d ^ oo, the order of a diagram has to be counted in a different way. Any nonlocal contribution becomes 
negligible, as is visible for example on the contributions of the third-order term of Eq. (^^ to the continued fraction 
Eq. (p?! ). Actually, the exact solution in d ^ oo is given by the diagrams depicted in Fig. H, as was mentioned in 
Ref. [43| . Therefore, the approximation (^^ becomes quite crude in infinite dimension. Note that our criterion for the 
Mott transition differs from that used in d oo, since we simply demand the closure of the gap at the Fermi level, 
without requiring the appearance of a new coherent peak in ^(k, cj) at a; = 0. Such a peak does not appear in the 
present approach. 

When T < u, we cannot calculate t^. analytically because we do not know which excitations close the gap, but 
Fig. H sketches a numerical evaluation (for d = 1) in the (T, t) plane of the line where the gap vanishes. The value of 
tc grows upon lowering T, and there is no Mott transition at zero temperature, in agreement with the exact result.u 
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FIG. 2. Diagrams yielding the exact r(ia') when d — > oo. 14uto(k,iaj) = V(k)/(1 — r(iaj)l/(k)) is itself dressed with these 
diagrams. 
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One has to keep in mind two important things about the insulator to metal transition described in this subsection. 
First, if we call Q the difference between the momentum of the highest negative energy (hole-like) excitation and 
that of the lowest positive energy (particle- like) excitation, then Q takes the value (tt, ...,7r) at high temperature, 
and goes to zero when T — > 0, but the gap always remains indirect. Secondly, even when the gap has been closed, 
the spectral weight at the Fermi energy still has several well separated peaks. A(k, uf) is never dominated at small 

by a single narrow quasiparticlc peak allowing to define a Fermi wavevector and a Fermi velocity. Thus, when 
extrapolating it beyond the transition, our solution does not describe a Fermi liquid. This is not surprising since we 
do not expect a perturbative method starting from the completely localized insulating state to be able to describe 
both phases satisfactorily. 
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FIG. 3. Crossover diagram of the one-dimensional half-filled Hubbard model as obtained from the third-order Green 
function. The estimated domain of validity of the approximation is the region under the dot-dashed line. 




C. Antiferromagnetic correlations 

Here again, a rigorous discussion of the magnetic behavior of the system should be based on the magnetic sus- 
ceptibility, which is a two-particle function. Nonetheless, the spectral weight should contain a valuable qualitative 
signature of increased antiferromagnetic correlations, namely a change in its periodicity. Indeed, if there were a true 
antiferromagnetic ordering of the spins, the doubling of the unit cell should lead to a periodicity of tt in reciprocal 
space, instead of 27r, in all the correlation functions. We expect signs of this reduced periodicity to appear, even in 
the disordered phase, when antiferromagnetic correlations become strong. Their effect actually show up at low T, as 
illustrated by the plot of A{k, uj) of Fig. ^, corresponding to point C of Fig. ^. For simplicity, we discuss the ID case 
in this subsection (k becomes k), but the conclusions apply to any dimension. We stress that we only observe the 
appearance of increased correlations but always remain in the disordered state: a true transition towards a long-range 
ordered state in d > 3 is not visible in our results, and Fig. ^ is by no means a phase diagram, but only a crossover 
diagram where we distinguished regions of different qualitative behavior of the one-particle spectral function. A{k, uj) 
has four delta peaks (a finite width 77 is added for clarity) given by dispersion relations ^^(fc), (i=l to 4, as shown in 
Fig.|). 

The spectral weight is an even function of k, and particle-hole symmetry at half-filling ensures that A(fc-|-7r, —u)) — 
A{k^uj). The minimum of L02{k) is at fc = at small t and high T (e.g. point A of Fig. |3|). But when T is lowered 
down to point C (Fig. ^), this minimum moves continuously from fc = towards k — 7r/2, and peak 2 loses weight 
for values of k much smaller than 7r/2. These changes reflect the AF short-range order that gradually builds up 
when T becomes smaller than the AF superexchange J = 2t^ /u of the equivalent t — J model. As suggested above, 
the approximate cell doubling in direct space translates into a nearly 7r-periodic dispersion for peak 2, although the 
27r-periodicity of its weight and of a;i(fc) reminds us that the state remains paramagnetic. This is the paramagnetic 
analog of shadow bands in the spectral weight of the cuprates. We chose to define an AF crossover line in Fig. || as 
the points where fc = ceases to be the minimum of the dispersion 012 (fc)- This crossover line is roughly parabolic 
(T ~ t^) at low t, which implies a crossover temperature proportional to J = 2t^ /u. Furthermore, in this regime, the 
width of band 2 is of order J whatever the value of i, supporting the above interpretation. 
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IV. COMPARISON WITH NUMERICAL RESULTS AND LIMITATIONS OF THE METHOD 



The purpose of this section is to try and find out where in parameter space - the (T, t) plane - the approximate 
Green function ( |2^ ) can be regarded as rehable, to improve it and compare it with Monte-Carlo results when possible. 
The entire section applies to the half-filled Hubbard model. 




FIG. 4. Spectral weight of the half-filled one-dimensional Hubbard model in the antiferromagnetically correlated regime 
(point C, t = 0.25m, T = 0.06ii = 0.24t). 



A. Reliability of the third-order solution 



It is not an easy task to define a domain of validity for the present approximation scheme for several reasons. 
First, the expansion parameter being the hopping amplitude t, any hopping process is considered on an equal footing 
regardless of whether it involves a change in the double occupancy or not. As a consequence, we expect the ratio t/T 
to play as important a role as t/u. Secondly the end result involves i in a complicated way, and not only through a 
power series. Nevertheless, since the atomic Green function is just 

-^-^ (30) 
u 

ILU 

iuj 

and gives spectral weight only at iu — ±u, we expect an approximation of the form ( p5| ) to be valid when the 6;'s are 
small with respect to u. When applied to Eq. (|27|), this criterion leads to the conditions: 

tV 1 fT 



2dt<u and - < TTl - ■ (31) 

\U J M \U J 

The first requirement is quite intuitive: it expresses that the bandwidth is smaller that the on-site interaction, which 
was the basic assumption anyway. The second one gives the low-temperature limitation of the method. In dimension 
d = 1, the region where the two conditions Eq. (^ij) are fulfilled lies under the dashed line of Fig. ^. When extrapolating 
the results outside this region, one cannot predict how fast the approximation may deteriorate. In particular, it is 
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line in 
results 



worth mentioning that the t — > oo (free-particle) limit is recovered properly. On the other hand, at small t the T ^ 
limit again gives a free-particle behavior (except for a singular behavior for wavevectors such that T^(k) = 0), which 
is obviously wrong. This allows us to define, besides the region of (almost) certain validity under the dashet 
Fig. |[ a hatched region of acknowledged failure. For the one-dimensional case, the theoretical^ and exactli^ 
describing spin-charge separation fall in this region of parameter space, which prevents us from making any meaningful 
comparison. However, a definite prediction of our work is that upon raising the temperature, there appears noticeable 
spectral weight near k = n (for ui < 0). This new feature of the spectral weight, visible as peak 3 of Fig. ^ and absent 
from the zero-temperature solutions, might correspond to the "question-mark" features in Fig. 1 of Ref. Thus, 
temperature seems to have a drastic effect on the distribution of spectral weight. 



B. Beyond third order 

Requirements such as Eq. (plf) are expected since we are dealing with a perturbative method. But beside these 
limitations, our solution (Eq. (P7|)) is plagued with a more serious problem, namely its discrete spectral weight. 
Indeed, the exact solution in the thermodynamic limit certainly has a continuous distribution of poles on the real 
axis. Such a continuous distribution is in general necessary to account for spin-charge separation, or a finite lifetime of 
one-particle excitations. Therefore, an approximation describing the spectral weight as a sum of four delta functions 
is necessarily a crude one. The reason why it is difficult to obtain a continuous spectral weight is the huge degeneracy 
of the unperturbed Hamiltonian. The starting point being a collection of independent atoms with the same two 
energy levels, it is not likely that a finite-order perturbation scheme can produce a macroscopic number of distinct 
approximate eigenvalues. We will discuss in section |^ a partially self-consistent approach leading to an extended 
spectral weight. But in the latter case, we no longer have the freedom to add higher-order terms in order to solve 
the causality problem. Facing these difficulties, we computed the fourth and fifth order of the Green function, in the 
hope that they would introduce many more poles into the two narrow Hubbard bands and give some hints of the 
lineshape towards which the spectral weight tends. Let us stress straight away that one is never certain to improve a 
perturbative result by adding higher-order terms. Actually, in most of the cases where perturbation theory is used, 
the series is asymptotic rather than convergent, and there is an optimal order - typically of the order of-the inverse of 
the expansion parameter - beyond which the approximation deteriorates quickly with each new termMS Computing 
the fourth and fifth order presented a substantial technical difficulty, since it involved the function G^^^*^ having a 
different expression for each of the 5! — 120 possible time orderings (one of the times can be set to zero). We overcame 
this difficulty by designing a symbolic manipulation program dedicated to this problem, taking advantage of the very 
systematic form of the expressions in the atomic limit. Let us mention that up to fifth order - and that seems to 
be general - the even orders lead to a modification of the partial numerators (a;) only, and the odd orders lead to 
a modification of the partial denominators (bi) only, within the continued fraction (p5[). Appendix |^ presents the 
diagrams and the result for r(k, iw) up to fourth order. r(k, iw) and the corresponding continued fraction have been 
computed up to fifth order included, but are too lengthy to be presented. 

At both orders the Green function has eight poles, among which only six have a significant residue. Thus, there 
appears only one more pole in each Hubbard band with respect to order t^. Furthermore, except for the region where 
the series certainly converges, the fifth order may give an answer quite different from that of the third or fourth order. 
In particular, the presence of diverging coefficients makes the self-energy (in the usual sense) go to zero very fast 
upon lowering the temperature, and the solution the n coin cides with the free-particle limit, completely missing the 



antiferromagnetic correlation effect described in Sect. IIIC Thus, despite the important effort invested in computing 
orders and , little progress has been achieved: the poles still appear as a sparse sampling of the spectral weight 
rather than a faithful fit, and the complexity of the solution has increased up to a point where it is hardly tractable. 

In order to test the analytic validity of our results, we computed the exact Green function up to fifth order for the 
exactly solvable two-site problem, and compared with the diagrammatic approach. The two approaches are proven 
consistent with each other and the test is summarized in Appendix In the precise case of two sites, the continued 
fraction has only four poles instead of eight, as it should. It is interesting to see that in that case too, the best 
approximation can be the third, fourth, or fifth order, depending on the parameters. Most importantly, the two-site 
problem confirms that dealing with a small value of i/u is not sufficient to ensure the accuracy of the approximation. 
As a matter of fact, we observed empirically that higher orders improve the solution for T > t, and deteriorate it 
for T <ti t. Explaining precisely why perturbation theory fails at low temperature appears difficult; this is probably 
related to the huge degeneracy of the system in the atomic limit. 
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FIG. 5. Spectral weight associated with points A,B,E,F,G and D of the crossover diagram (Fig. g). For each point, we give 
the spectral weight obtained from cpder (top left), (top right), and (bottom left). A finite width 77 = 0.02 was added 
for clarity. The Monte-Carlo resultal3 (smooth curves) were used to establish the density plot (bottom right). Note that: 
Point A lies within the insulating paramagnetic region, as can be deduced from the presence of a gap and the monotonous 
dispersion of the bottom of the upper band (see text for details). For these values of the parameters, order seems the best 
approximation. 

Point B is in the antiferromagnetically correlated region: The minimum in the dispersion of the bottom of the band has shifted 
to fc > 0. Orders and are best suited there. 

Point E is in the insulating paramagnetic region and is best described by the fourth-order result. 

Point F is just at the edge of the antiferromagnetically correlated region. The third-order approximation is the best one. 
Point G is in the metallic region. Only order gives the closure of the gap for this point. 
Point D is well within the metallic region. 
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C. Comparison with Monte-Carlo data 



We now present Monte-Carlo data supporting and illustrating the conclusions claimed so far, regarding the utrysical 
behavior of the model, as well as the accuracy of our approximation at various orders in t. The simulationafj were 
done for a twenty-site one-dimensional lattice for reasons of computing power a*ijd time. The spectral weight was 
deduced from imaginary time Green functions by the Maximum Entropy MethodB3 The first interesting region is the 
crossover between the ordinary paramagnetic insulator and the strongly AF correlated insulator. Fig. |^ shows the 
spectral weight for the points A {t = 0.25m, T = 0.25) and B {t = 0.25u,T = 0.125) on each side of the crossover. 
Third-, fourth-, and fifth-order results are displayed. Of course, A and B being ve ry clo se to each other, they have 
largely similar spectral weights, but their qualitative difference defined in subsection III C is supported by the density 
plot. The fact that orders and t* are equally well suited for these two points, located respectively on the lines 
T = t and T = t/2, is consistent with our argument that Order deteriorates compared to Order precisely in that 
region. Let us mention that in one dimension, the Monte-Carlo data show that the antiferromagnetically correlated 
region is tiny, and that a spectral function reminiscent of spin-charge separation (with an important transfer of weight 
from low to high energy for k < 7r/2, a; > and k > 7r/2, oj < 0) appears when lowering T a little further (actually 
when entering the shaded area of Fig. |3|). However no such thing is expected in two dimensions, for which the 
antiferromagnetically correlated region extends without spin-charge separation down to T = where true long-range 
order takes place. 

The upper border of the paramagnetic insulator region, defined by the closure of the gap, can be computed more 
precisely than in Fig. |^ with the help of the fifth order. It leads to a lower value of tc, saturating around tc ~ 0.47u at 
high temperature. Order has been used when establishing the improved crossover diagram of Fig. a The sequence 
of points E, G, D, taken along the line T = 1.25t for t = OAu, O.bu, 0.8u, and also depicted in Fig. |5, shows that the 
gap actually closes for a smaller value of t than predicted by order 3. When lowering the temperature below t, the 
fifth-order result quickly deteriorates: it does not yield the antiferromagnetic crossover, and predicts the closure of the 
gap for smaller and smaller t^, with the wrong conclusion that tc when T — > 0. Thus, in Fig. 0, we use the third- 
order result to obtain the antiferromagnetic crossover line at low temperatures, and we switch from order to order 

at the crossing point between the insulator to metal line (at high T) and the antiferromagnetic crossover line (at 
low T). Therefore, the position of the line separating the paramagnetic metal (gapless) from an antiferromagnetically 
correlated insulator (gapped) becomes conjectural, since fifth order fails to describe it, and third order no longer 
brings it to the right "triple" point when lowering t. Anyway, the density plot for point F {t = OAu, T = 0.35) (Fig. ||) 
clearly shows the tendency to open a gap upon lowering the temperature. 




FIG. 6. Spectral weight deduced from our approximate Green function at order for Point H2 of Fig. ^ (t = Pt25u, 
T = 0.025ii) in dimension d = 2. This is to be compared with the Monte-Carlo data of Ref. As in dimension d — 1,e3 our 
result agree quite well with the numerical computations despite the low value of the temperature. 
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FIG. 7. Improved crossover diagram of the half-filled Hubbard model in dimension d = 1 (top) and d — 2 (bottom). The 
dot-dashed line reminds the estimated validity region in one dimension. The antiferromagnetic crossover was calculated with 
the order result, and the Mott transition line with the help of order . The limit between the metallic and the insulating 
antiferromagnetic regions is no longer clearly defined. 



Finally, Fig. || shows a comparison of the solution and the Monte-Carlo data of Ref. for point H2 in the 
two-dimensional case. The agreement is excellent, even though the parameters fall in a region where we do not expect 
our approximation to be good. Such a good agxeement with Monte-Carlo data outside of the expected region of 
validity has also been observed in one dimension. Ej This probably means that the criteria are to stringent. One 
may instead apply the following a posteriori criterion: The approximate spectral distribution is reliable if a large 
fraction of the weight (say, 90%) falls within ±t of the original Hubbard bands at w = ±u (this is the case for Fig. ||, 
as well as for Fig. 3 of Ref. ^). This softer criterion of reliability probes the overall effect of hopping on the spectral 
weight, without a detailed analysis of partial numerators and denominators, whose effect on the spectral weight is 
often far from obvious. 

To summarize, the strong-coupling expansion has two serious limitations: (i) one cannot expect to describe faithfully 
the lineshape of the spectral function and (ii) it is impossible to obtain an accurate result for any temperature much 
lower than t. We add that going beyond fifth order within the systematic approach described in this paper would 
be difficult and unpractical, since the (intractable) result would be relevant only deep in the insulating region. On 
the other hand, Monte-Carlo calculations agree well with our results as far as the overall distribution of the weight is 
concerned, and confirm the qualitative conclusions of section III, summarized in Fig. ^ Thus, we are confident that 
the improved crossover diagrams of Fig. |7| are reliable. 



D. Double occupancy 



The double occupancy {n^-ni) is a local static quantity essential to the comprehensioB-|Of the Hubbard model, which 
gives in particular the average potential energy U {n-^ni) . Except in infinite dimension,Ell (n-|-n|) has to be determined 
with the help of Monte-Carlo simulationsHH It is possible to deduce it from the one-particle Green function through 
the exact relation : 

-i^ *^)S.(k, iLo) e«^ = U{n^ni), (32) 
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where So-(k, iuj) is the self-energy. There are two possible ways for us to use the identity (32). On the one hand, from 
our knowledge of the thermodynamic potential fl (cf. Sect. [V| below), we may extract the power series in t (thereafter 
denoted St) of (n^ni): [3st — —dfl/dU. On the other hand we can insert the continued fraction ( |27| ) into Eq. ( |3^ ) 
and compute (n|7i|) numerically This subsection shows that the second option (thereafter denoted sj) gives much 
better results, which confirms that the continued fraction representation is a controlled way of resumming diagrams. 
At half-filling, the power series for the double occupancy is: 



1 / ^'e^(-l-f e^) ^e^ -1 + e^ . 

^'~l+el^^ [ (l+e^)3 " (1 + e/3)2 + 2(1 + e/3^ ' 



2 



f 9{-5 + d){-l+ c^) 3(4-H7e^-3de^-H4e^^) , 2 (l + 7rf)e^(-l-H e^) 
+ V 32(1 + e/3) ^ 16(1+ e/3)2 16(1+ ^3)3 

^3 3(„l + 3rf)e/3(l_4e/^+e^/3) + 7d) e^(-l + e/^)(l - lOc^ + c^^^) ^ 

^ 8il + W 16(l+e/3)5 (33) 

where we have set u = 1. Fig. ^ sketc hej . % and sj as functions of t along the line T — t/6 in the two-dimensional 
case, as well as some Monte-Carlo dataHa Whereas St quickly increases beyond the maximum physically acceptable 
value of 1/4, sj interpolates quite well between the low t regime and the high t regime where the form of the continued 
fraction ensures that the free-particle limit (n-^ni) — > 1/4 is properly recovered. The effect of temperature on double 
occupancy is summarized in Fig. |^, showing sj as a function of T for several values of t, as well as Monte-Carlo results 
from Refs. p2 and ES. 
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FIG. 8. Double occupancy (n^ni) as a function of the hopping r^rameter along the line T = t/Q in dimension d = 2. 
Comparison between St (dashed), sj (solid), and Monte-Carlo resuItsH (diamonds, bigger than the uncertainty). Both st and 
sj are calculated at order t^. 
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FIG. 9. Double occupancy (nin^), computed from sj at order t*, for d = 2 as a function of temperaturefJor various values 
oi t (t = 0.5m, t — 0.33u, t — 0.2u from upper to lower curve). The diamonds represent Monte-Carlo resultsO 

sj always presents a minimum at some temperature, which may be more or less sharp. Such a minimum is 
obtained in infinite dime*i6ion,Lj where it is explained by a Pomeranchuk effect. A recent application of a two-particle 
self-consistent approached shows that this minimum is still present in dimension three, but disappears in dimension 
two: In low dimensions the self-energy acquires a momentum dependence that makes it more sensitive to magnetic 
correlations, thus suppressing the Pomeranchuk effect. In any case, we believe that this minimum in sj is an artifact 
of our approximation, not an entropy effect. Indeed, when lowering T well below that minimum, we enter the regime 
where the self-energy artificially goes to zero with temperature, and consequently the double occupancy goes to 1/4. 
We have stressed several times that this limit is not well rendered by our solution, and it is confirmed by Monte-Carlo 
data which show little T dependence at low T. Anyway, as long as we stay in the vicinity of this artificial minimum, 
our theoretical predictions differ from the Monte-Carlo values by at most 15% for t — 0.5u{U — 4t), and the agreement 
improves quickly at high temperatures, even for values of U in the intermediate coupling regime, as in Figs. I and |. 
It is expected that the agreement will be better in higher dimension. Moreover, at low temperature, the agreement is 
better for smaller values of t/u, as expected. 



V. INFINITE RESUMMATIONS 



Up to now, we have followed a systematic approach, which consists in building the series in powers of i of a given 
quantity. This has led us to interesting results, but has also shown its limits: (i) it does not yield continuous spectral 
functions and (ii) we have pushed the calculations to their maximum acceptable complexity level. At this point, rather 
than keeping all diagrams of a given order, we look for infinite series of diagrams that might embody an interesting 
physical feature of the model. 
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FIG. 10. The one-loop "rosary" diagram. The simple hopping term is replaced by an RPA-like propagator. 
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FIG. 11. Spectral weight obtained with the auxiliary self-energy of Fig. with hopping term t — 0.5m, and momentum 
k = 0.57r (full curve) and k — OAn (dashed curve). The spectral weight is not everywhere positive within this approximation. 
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A. RPA-like propagator 



The most natural possibility is to resum all the 
gator t^(k) by 



'rosary" diagrams, that is, to replace the auxiliary fermion propa- 



F(k) 



1 - V(k)G{iuj) 



(34) 



wherever it appears in the diagrams. This prescription corr espo nds to an alternate way of separating free and 
interacting parts of the action, as mentioned at the end of Sect. II B . With this new propagator, the electron, instead 
of simply hopping between two neighboring sites, can travel an arbitrary sequence of one hop (of amplitude VCk) 
followed by a period of rest on some site (of amplitude G{iuj)). The simplest series of this kind is shown on Fig. |lC 
The result in d = 1 is: 



lU! 



3^2 



-1 



V 



I-, _ f 2tiuj \ 



(35) 



Let us comment on T{iuj), which coincides with the entire Green function when k = ■n/2. Our main objective has 
been fulfilled, in the sense that the spectral weight A{t:/2^uj) is now a continuous distribution in the band 



- t + + <UJ <t + 



(36) 



and is symmetric. On the other hand, y4(7r/2, lo) is not everywhere positive. Indeed, a detailed analytic study of Tiiuj) 
shows that it has a positive spectral weight everywhere except at a; = ±u. At those points, the first term of Eq. (pq) 
brings a delta function of weight 1/2, and the second term a delta function of weight —3/2, resulting in a negative 
peak in the spectral weight. The prefactor 3 in front of the second term has been thoroughly checked, and we have 
achieved great confidence that Eq. (35) is correct. When k ^ 7r/2, the spectral weight has to be derived from the 
imaginary part of 



Q{k,iuj) 



1 



r(iw)-i + 2icos(fc) 



(37) 



when iuj —t Lo + irj. Like for k — n/2, A(k,Lu) is nonzero within the bands, and negative ai lu — ±u. In addition, 
A{k,U!) has two delta peaks near the edges of, and outside the bands - where the Green function has no imaginary 
part. Fig. |ll| shows A{0.5n,uj) and A{OAtt,uj) for t = 0.5m. 

Thus, the one-loop diagram of leads to a normalized, but nonpositive spectral weight. We have also lost the 
interesting physical effects of order three: the result is temperature-independent and always has a gap at the Fermi 
level. We did not find any acceptable way to cure this negative-weight problem. One could suggest to reduce the 
normalization of the one-loop diagram of Fig. pXl, but introducing such an arbitrary factor is by no means justified. 



k = n 




(X u) 



FIG. 12. Spectral weight of the one-loop self-consistent solution, for t — tc — 0.39u. The wavevector goes from fc = (top) 
to A: = TT (bottom). The spectral function is clearly not normalized within this approximation. 
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B. One-loop self-consistent approximation 



Another natural attempt involving an infinite subset of diagrams is to compute tlie one-loop self-consistent solution: 
We keep the first two diagrams of Fig. ^, and discard all vertices beyond G^^'^. The self-consistent solution f (io;) obeys 
the following equation: 

f.M = G,M--i^ ^ Gi,i;^ -. (38) 

, — . 1 — K ki i iwi 



ai ,ki ,iuj 



The sum in Eq. ( p8| ) can be computed exactly, given the following property of G^^'^ at half filling and zero magnetic 
field: 

E^TTr —3f3u^S{iuj — iui) 
G^,<Ti;<T,^l(*^>*^i;*^,*^^l) = 7-—Z Z-2 +t:{lLU,lUl), (39) 
((iw)2-u2) 

where £ is an even function of iuji. Since we expect T„{iu!) to be antisymmetric in icj, the self-consistent propagator 
in Eq. ( ^8[ ) is antisymmetric when icoi — > —iuji and ki — > (tt, ...,7r) — ki. Therefore, £ does not contribute and 

Eq. (^^ can easily be solved numerically. Among the various solutions for the spectral weight, we retain the positive 
one having a compact support at very small t, and follow its evolution when increasing t. The result is shown in 
Fig. [ij for t = 0.25w, and is obviously not normalized. The quick suppression of spectral weight close to a; = ±u 
when k differs slightly from 7r/2 is a surprising feature too. However, we have recovered the closure of the gap for 
t ~ 0.39m, a value close to the prediction of Order at high temperature. 

The two examples just described show that it is not easy to include self-consistency in strong-coupling perturbation 
theory. Lack of positivity or normalization appears in the simplest attempts, and there is no clue on how to solve 
this problem. However, it seems the only way of obtaining a continuous spectral function, and further developments 
should probably include a dose of self-consistency. 



VI. DOPING 



Studying the Hubbard model at arbitrary filling is important for several reasons. First of all, the charge gap of 
the half-filled system disappears upon doping, and in dimension d — \, bosonization shows that this implies dramatic 
changes in the spectral weight. Furthermore, in dimension d = 2, a slight doping is directly relevant to the high-Tc 
superconductors. Finally, the metal-insulator transition can be induced by doping rather than by interactions. We 
address the latter aspect in this section. 

On technical grounds, working with arbitrary chemical potential is extremely difficult, since the complete expression 
of G^^'^ given in appendix ^ is already hardly tractable, and using it to compute diagrams would be even worse. Even 
if it is obviously preferable to include the full atomic Green function in the unperturbed part, one could treat the shift 
in chemical potential = n — u as a. perturbation, as is done in appendix ^ However, the spatial integrals would 
no longer cancel several low-order diagrams involving G^^^*^ and G^^*^. Hence the calculation was carried out with a 
given /i, without expanding it as — u + t^, with the help of a suitable generalization of our special purpose symbolic 
manipulation program. We were able to compute F up to order in dimension d = 1, and to build the corresponding 
continued fraction. The partial numerators and denominators of the latter are given in Eqs. (|4^,^) below. 

1 

-4u[ - 2/?(l - 2i/)2(i/2 - V2)t^ + (-1 + v)vu + 0^{l - 2vf x {-v + Av"^ + V2 - 81^1^2 + 'Avl)t^u\ 
^^2 I Af3cos{k)i-l + 2i^){iy-iy2)t^ 

-4(-2-z/ + ;/2)y2 8/3cos{k){-l + 2i^f{iy-i^2)tu^ 

9 21{-l^v)v ^ ' 
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bi = -2tcos(fc) +to- (-1 + 2;/)[4/3(-i/^ + 1/3)^^ + u + 2l3'^{~iy + Ai^'^ + V2- ^vv2 ^ -ivDt^ u] 



ho = 



b3 = 
where 



to + (-1 + 2iy)u + J^^—L^^\iy^ + v\-2 + AI3u) + i/3(2 - /3(5 + 8^/2)^) 
(-l + z/)i/ L 

+ 4/3i/|m + 1^2(2 + Pu)) + z/2(/3u + ^pj^lu + 1/2(2 + 9/3w))} 

WcoB{k)t^ + 4/3u) + z/2(-l + Piy2u) + v\l + l3u + AI3vlu 

(—1 + v)vu L 

+2z/2(l + 4/3u)) - 2t/i/2(l + /3(u + 2v2u)) - 2v^{-l + 2/3(m + 2i/2u))} 

3^0 + ^ -2;/^ + 2cos(fc)^ r^^^ ^ ^ ^ ^ _ ^^^3^ ^ ^ _ ^ ^ 
3 9(— 1 + vji' L J 

to + ^'^t^''^'' + I - 4/3Z/2M + 1/2(9 - 16/3(1 + V2)u) + IG/Ji/^w + i/(-9 + 4/3(m + 4z/2u))| , (42) 

3 9( — 1 + 1/):/ L J 



2/3*0 1 p/3(to+") 

+ (43) 



1 _l_ e2/3*o _|_ 2 e'3(to+") 
is the average number of electrons of a given spin per site in the atomic hmit, and 

" 1 + c2/3to +2e/5(*o+'') '-^^^ 

the double occupancy in the atomic limit. Although complicated, Eqs (|4l| , |4^ ) are plausible for several reasons. 
First, our symbolic manipulation program was checked on the thermodynamic potential, and our result agree with 
Rcfs. 32 3^, indicating again that the result of Ref. |3l| is incorrect. Secondly, the half-filling case is properly recovered 



from Eqs. (|4T]j4|) 

At high-enough temperature, the spectral weight derived from Eqs. ( |4l|j42| ) evolves smoothly with the chemical 
potential. For example, for < io ^ the Fermi level shifts slightly towards the positive energy peaks, whereas the 
weight of the latter increases slightly, with respect to the half-filled case. A similar conclusion applies to the density 
of states : Increasing /i shifts the Fermi level towards the upper Hubbard band, and gives the latter more weight, the 
overall effect being an increase in the occupation number. In other words, except for a slight redistribution of weight, 
the behavior of the system resembles that of a band insulator. 

However, this smooth picture collapses when the temperature is too low, the spectral weight abruptly becoming 
negative at various energies, which means that the a{s of Eq. ^ are no longer all positive. This is yet another 
manifestation of the limitation of our method at low temperature where, according to Rcfs. 24 2^, there should be a 
quick and massive redistribution of spectral weight between the Hubbard bands when varying /i. The breakdown of 
our solution, which occurs around T ~ is concomitant with a lack of monotonicity in the relation between chemical 
potential and filling. This relation is implicitly defined by the thermodynamical potential through the relation 

and becomes ambiguous as soon as the thermodynamical potential is only known approximately. From this expression 
a truncated power series in t for n (denoted n'--^\fj,) at order t^) follows directly. Alternately, one may reverse this 
relation and express /i as a truncated power series in t, which we write /x^''-' (n) . The two expansions n'-'"' (/x) and /x'^''-' (n) 
may lead to different physical conclusions, which would be unacceptable. We have verified that, when lowering the 
temperature, both functions n'^'^^{p) and fi^^\n), calculated at Order i^, cease to be monotonous nearly at the same 
value of T/t, which is also the temperature at which the spectral weight becomes unphysical. 

Nonetheless, and similarly to what we did in the half-filled case, we can use the penetration of the Fermi level 
into one of the Hubbard bands as a qualitative criterion for the metal-insulator crossover. This penetration should 
translate into an important increase of occupied states (likely to be delocalized) in either Hubbard band, and therefore 
an increase in the conductivity. Fig. |l^ shows, for several temperatures, the line in the {t,n — l)-plane where the 
Fermi level enters the upper Hubbard band. The numbers associated with each curve are estimations of the accuracy 
at the points T — t, beyond which the solution breaks down. These numbers are the relative difference between the 
results for n(/x) obtained in two independent ways : from the Green function (through the density of states) and from 
the thermodynamical potential. 
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Note from Fig. |l3| that, in the atomic hmit t ^ 0, the transition occurs near n — 1 = ^. Indeed, in the atomic 
Umit, it is a simple matter to show that the chemical potential vanishes at T = at this filling, and finite temperature 
corrections are exponentially small as long as T < u. The chemical potential then touches one of the (infinitely 
narrow) atomic bands (at tu = and cu = 2u) and any infinitesimal value of t makes the system conduct. 




FIG. 13. Location of the metal-insulator crossover at various temperatures. The numbers are estimations of the relative 
error on n — 1 at the points t = T, given that the solution quickly breaks down for t > T. 

Thus, in the high-temperature regime where our method is reliable, the Mott-Hubbard insulator shares som of the 
features of a band insulator : the gap is nearly rigid, and a crossover towards a good conductivity is obtained when 
the population of carriers ceases to be exponentially small. In contrast to a band insulator, there is, however, some 
transfer of spectral weight towards the band closest to the Fermi level, when /z is changed from its half-filled value. 

For T < t, the expected drastic changes in the spectral weight are not properly described by our approximate 
solution, which instead becomes unphysical. 



VII. CONCLUSION 

In this work, we have developed the ideas of our earlier paper,S and pushed the computations to higher order. 
The initial causality problem, a conspicuous difference between strong- and weak-coupling perturbation theories, 
was solved by the continued fraction representation. Among our main results are the crossover diagrams of Fig. |^. 
They demonstrate the possibility of describing the three qualitatively different regions of parameter space of the 
half-filled Hubbard model by a single approximate analytic expression for the Green function. With respect to 
Ref. these crossover diagrams were refined by higher-order terms. Moreover, the crossover diagrams are based on 
spectral weights that we have shown to be consistent with Monte-Carlo simulations. The difference between these 
spectral weights and their zero-temperature counterparts obtained from other methods, demonstrates an important 
temperature dependence of A(\i,Lu). 

We have pointed out the limitations of the simple strong-coupling perturbation series. Firstly, our solution breaks 
down quickly upon lowering the temperature just after the magnetic correlations have started to show up. This 
explains why the spatial dimension d plays a minor role in our solution. Actually, as long as magnetic correlations 
play no role, i.e., at high-enough temperature, all dimensions are known to show similar behaviors as far as the 
metal-insulator crossover is concerned (we include dimension d = 1, for which the gap predicted by Bethe Ansatz is 
exponentially small in a finite range of U/t). But our method fails at lower temperature, where the system evolves 
towards spin-charge separation {d = 1), or a true antiferromagnetic transition at zero {d = 2) or finite (d > 3) 
temperature. Secondly, the systematic strong-coupling expansion can only provide a finite number of poles in the 
continued fraction and the correponding approximate Green function cannot have continuous spectral weight. This 
does not prevent meaningful comparisons with Monte-Carlo data, since the overall spectral weight distribution can 
be assessed. This "discieteness problem" might be solved with the help of a suitable termination function within 
the continued fraction,Efl but this is impossible without any prior knowledge about the Green function. The simplest 
ways of obtaining a continuous spectral weight from infinite subsets of diagrams lead to nonpositive or unnormalized 
functions. 

Another important difficulty of any expansion about the atomic limit is the absence of Wick theorem for the 
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atomic limit itself. The auxiliary fermions allowed us to organize the expansion in the best possible way, and to 
reduce tremendously the number of diagrams at a given order. Nevertheless, the exact G^'^'s are still necessary, 
and the result of appendix ^ for G^^"^ suggests that they are quite difficult to compute for i? > 3. Introducing 
a physically motivated approximation scheme for the atomic vertices would allow further progress, particularly on 
the two-particle correlation functions. The full two-particle correlation function in the atomic limit {U — > cxi) is a 
noteworthy byproduct of our work. 

The major challenge faced by the strong-coupling perturbation theory is the low-temperature barrier. We pointed 
out that the exponentially large degeneracy of the atomic ground state is likely to be the source of the problem. 
One possible solution to this difficulty is to select a ground state more likely to connect with the low-temperature 
phases and to organize the perturbation series around that ground state, which would imply a certain amount of 
self-consistency. Work along these lines is in progress. 

While this paper focused on properties derived from the one-particle Green function, two-particle Green functions 
are also accessible within the method presented here, but their systematic computation is more involved at order 
(it requires the atomic four-particle function G^^). For this reason, we defer discussion of two-particle correlations to 
a future publication. 



VIII. ACKNOWLEDGEMENTS 



We thank C. Bourbonnais and N. Dupuis for many useful discussions. We are grateful to H. Touchette, S. Moukouri, 
L. Chen, and especially D. Poulin and S. Allen for sharing their numerical results. Monte Carlo simulations were 
performed in part on an IBM SP2 at the Centre d' applications du calcul parallele de I'Universite de Sherbrooke. This 
work was partially supported by NSERC (Canada), by FCAR (Quebec) and by a scholarship from MESR (France) 
to S.P. 



APPENDIX A: DIAGRAMMATIC RULES 



This appen dix is devoted to deriving and illustrating the diagrammatic theory valid for the auxiliary field introduced 
in Sect. II B. For definiteness, we suppose that the site index describes a d-dimensional hypercubic lattice, and we 
only consider Hamiltonians where the hi's do not explicitly depend on i, so that the interaction terms are translation 
invariant, in addition to being local in space. 

The first step consists in expanding the exponential of the interaction terms of Eq. (p^: 

z= /[d^'^ci^]e-^«[^^^]5:LLK:^,^J^^^] . (Al) 

For a given P, and taking into account the factor (— )^/P!, we have a sum of terms of the following form: 

(-^"^ -Sf^\i,\iA-SflW.^\. (A2) 



Cx\..Gh\ 

where Ri,..,Rp are P integers, and Ci,..,Gh the multiplicities of the different values that occur in the sequence 
.., Rp. Suppose we are interested in the following correlation function: 

Vfo^o ^{-f"UaO-^a»rbO -^o)- (A3) 



The contribution of a given power P to (A3) is the sum, for all possible sequences .., Rp, of 



Z 

p=l..P 
r-=l..R„ 



times a overall factor 



C,l.CHKRil)^..iRpir 



(A5) 
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which takes into account the minus signs coming with S^^^. Again, let us stress that the indices in the sum do not run 
over all their possible values, but rather are restricted to have a common site index if they refer to the same vertex 
(this restriction is encoded in the primed sum). Since the action is Gaussian, Wick's theorem is valid, and one can 
compute the average in expression (A4) as a sum over all possible permutations oi Rq + Ri + ... + Rp elements: 

'd {a^ } Gauss \ / Gauss 

p=l..p 

r=l..fip 

where (— stands for the signature of the permutation, and {i' a' "^^(^i,!)) Gauss ~ ~^a^-a{b^)- ^ term in the sum over 
the permutations can be represented easily by a Feynman diagram, according to the following rules: 

1. Draw P polygons (vertices) having 2Ri, ..,2Rp apices (internal points) respectively, and 2i?o isolated points 
(external points). 

2. To half of the external points and half of the apices of each vertex, attach an outgoing arrow (a "bra", cor- 
responding to a i/i). To the rest of the points and apices, attach an incoming arrow (a "ket", corresponding 
to a ip*). 

3. Label the bras in a conventional order (first the external ones, then all the bras of each vertex consecutively), 
and label the kets with corresponding primed integers. 

4. To each point assign a latin index, the indices of the internal points sharing the same value of the spatial 
component for each vertex. 



5. Finish drawing the arrows according to the permutation i!}. To the lines assign a propagator Vab (bra index 

R c 

first), and to the vertices assign a G^j^p^^^p^ (ket indices first). 

6. Integrate over all the internal variables. 



A few remarks are in order. First, the ratio Zq^uss/Z allows a restriction to connected diagrams onlyEJ Secondly, 
it can easily be seen that the exchange of two equivalent vertices does not change the value of a diagram, nor does 
the exchange of two internal points of the same type (bra or ket) on the same vertex. This allows a restriction to 
topologically different diagrams, each of which with a proper multiplicity factor. Thus, to span all the permutations, 
it is sufficient to consider only the topologically distinct diagrams, and to assign them a factor 

(^-)Ri+--+Rp+-^ X multiplicity 

Ci!..Cff!(i?i!)2..(i?p!)2 • ( ^ 

Due to the presence of various types of vertices (n-body interaction, n = 1,2,...), the counting of the sign and 
symmetry factor cannot be simplified as in the usual two-body interaction case. A third remark is that in general Ti" 
and Ti.^ conserve spin, and so spin is conserved along the lines, and the total spin entering a vertex is the same as the 
total spin leaving it. 



Qj, k^, i(0^ 




iaH-i(02-ia)2 



FIG. 14. Diagram of order \Vf contributing to the Green function of the auxiliary field Vc7(k, io;). 
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Oj, kj, icoj 




FIG. 15. Diagram of order |V^j* contributing to the thermodynamic potential. 



After a Fourier transform, space and time translation invariance of the on-site Hamiltonian h result in the conserva- 
tion of total momentum and total Matsubara frequency at each vertex. The rules remain the same, except that one can 
a priori impose the above conservation laws along the lines and at each vertex. There remain Ri + + Rp — Rq — P +\ 
independent internal variables entering sums of the following type: 



— y 



(A8) 



where L is the size of the lattice, k a vector of the reciprocal lattice, and ik^ = i{2n + 1)ttT a fermionic Matsubara 
frequency. Finally, if the on-site Hamiltonian h is sufficiently simple as to allow for an exact computation of the vertices 
(the G^'^'s), we can obtain, by usual diagrammatic perturbation theory, any correlation functi on o f the auxiliary field. 
We still have to derive from them the electron correlation functions: this is done in Sect. II C, Fig. shows an 



example of a diagram of order 6 contributing to Vab = ^{'^^a^'D- The algebraic expression corresponding to Fig.n4 is 



S 



J2 l^(ki)t/(k2)F(k3)^^(k + ki-k2 



xG 



lie 



IIIc 

(o'4),0-2,0'3;CT,(T3,(T: 



(iu) + iui — iuj2, iuj2, iuJs; ii^, 1^3, iwi), 



(A9) 



where S is the symmetry and sign factor. The value of can always be inferred from that of the other spins. In the 
present case, i?i -I- i?2 = 2 -|- 3 = 5, so the global sign is opposite to that of the permutation between the unprimed 
and the primed indices: —1. The symmetry factor is determined by the numbers between parentheses: the first one is 
the order in which the lines were drawn, and the second one the number of equivalent possibilities of drawing them. 
The result is 



S 



2x3x3x2xlxl_ 1 
(1!)(1!)(2!)2(3!)2 



The rules for the thermodynamic potential are nearly the same. From Eq. (0) one has 



VL = n 







riog / [dV'^dV']e('^l^"'l'^)(Je<'^l^>+<''l'''>^ 



/o- 



(AlO) 



(All) 



We know fJo exactly, and the second term is obtained by the sum of all connected diagrams with no external points, 
times — T. If we consider the thermodynamic potential per site, the factor becomes —l/{l3L'^). When computing 
the diagrams in momentum space, the factor l/{f]L'^) is already taken into account by the definition of the various 
Fourier transforms. Fig. |l5| shows a diagram contributing to the thermodynamic potential at order 4. The algebraic 
expression corresponding to FigJlq is 



J2 V^(ki)nk2)V^(k3)F(ki-f k2-k3) 
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^<^"%2;<T3,(<T4)(*'^i'*^2;«W3,zu;i +iuj2 - ii^s) 

(<T4);cri,<T2 (^'^Sj + *'*-'2 i 1 ^^^2 ) • (-^12) 



Here, i?i + i?2 = 2 + 2 = 4, so the global sign is opposite (see Eq. (All) and the following paragraph) to that of the 



permutation between the unprimed and the primed indices: —1. The symmetry factor is 



2xlx2xl_l 

(2!)(2!)2(2!)2 ^^8- ^ 



APPENDIX B: ATOMIC CORRELATION FUNCTIONS 

This appendix gives the one-particle and two-particle correlation functions of the atomic limit of the Hubbard 
model. The chemical potential is /i and there is a uniform magnetic field h. 

Hat = C/c|cJc4C| - ;u(c|c| + cjcj - /i(c|c| - c'^ci). (Bl) 



1. Green function 

The Green function 

GAri,T2) = -(r,c,(Ti)4(T2)) (B2) 
is straightforward to obtain. Introducing the mean occupation for each spin 

"± = 1 + Q{l^+h)l3 Q{,ii-h)l3 e(2M-C/)/3 ' ^^^^ 

the Fourier transform of the Green function reads: 

. + . ^ (- = ±)- (B4) 

io; + /i + ±/i luj + ^ + ±n — U 

2. Two-particle function 

We start from the definition of G^^ in imaginary time (by enclosing 174 between parentheses, we mean that its value 
is to be inferred from that of the other three spins): 

G'"^2,(<T4)<T3(n:T2;T4,r3) = (T^C^i(ri)c^2(T2)43(T3)c|^^)(T4)^ (B5) 

and compute it explicitly for any time ordering, using e^^^^^^^"^ as evolution operator and inserting complete sets 
of states as necessary. We then calculate its Fourier transform, defined by the following equation: 

f dn...dn G!,\,,.(.,),3(ri, r^; r4, r,) e-i-i+-^— 3x3-^^4x4 ^ (B6) 
Jo 

f3S{iuJi +iuj2 - ii^3 ~ it^A) G]J^^^^(^^^)^^{iuji,iuj2; {iuj4),i(^3)- (B7) 

Again, ilu^ = iuji + iuj2 — ito^. After introducing the following short-hand notation 

2 = 1+ e(^+'^"=' + e(^-''''3 + e'^^^-'^)'^, (B8) 
xf = iujj + fi±h, (B9) 



xf = iujj +ti±h-U, (BIO) 
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the connected correlation functions read: 



and 



G 



lie 

iT.iT 



{iujiiuj2, (10^4)1^3) = 

n+ + n_ — 1 
iuji + iuj2 + 2/i — [/ 



H h H h ' 



n+ — n_ 


M - 












2/1 Va^r 




( 2:4 






(3U^S{iuj2 - 
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n_ - 1 
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X2 
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-^2 ■^S 


1 '^4 



(Bll) 

(B12a) 
(B12b) 

(B12c) 

(B12d) 

(B12e) 

(B12f) 



All values of G^J^^^ ^^^-^^^{11^111^2, {1^4)11^3) can be deduced from Eqs. (Bll) to (B121) above, the antisymmetry with 
respect to the exchange of two incoming or outgoing variables, and the possibility to exchange t and J, by doing 
h ^h. 

Apart from the r egular terms of Lines (B12c) and (|B12i[) , G^^'^ contains only singular or possibly singular terms. 
Indeed, Line ( B12d|) h as 6{iu!2 — iuis) in factor. Line ( p312b| ) becomes proportional to Siiui + 1002) in the half-filling 
limit, and Line (B12c) becomes proportional to 5{iL0i — ilu^) in the zero magnetic field limit. These special cases have 
to be properly taken into account while computing G^^*^ directly at half-filling and zero magnetic field. 




FIG. 17. Spectral function of the atomic model with chemical potential fi — 1.8u at temperature T — 0.5m. Comparison 
between the exact result (solid), the raw approximation (dashed), and the continued fraction (dot-dashed). Whereas the dashed 
curve is neither normalized, nor even positive, the difference between the other two curves is hardly noticeable. 
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FIG. 18. Diagram contributing to F at fourth order (top row) and fifth order (bottom two rows). 



APPENDIX C: A PRACTICAL TEST OF THE METHOD 



In this section, we test the diagrammatic theory and the continued fraction representation on an exactly solvable 
simple model. We consider the atomic limit itself, with a slight doping given by a chemical potential ^ — u + Iq. On 
the one hand, we know the exact Green function from Eq. (B4): 

1 — T? Tl 

GexM= ■ ^ +-■ -T^ (CI) 



with 



g(«+to)/3 _|_ g2to/3 



(C2) 



On the other hand, if wc consider the chemical potential shift as a zero-range hopping perturbation, we can compute 
F to first order in powers of to with the help of the diagrams depicted in Fig. and obtain 

^ i^"^) ^ J-^ 2+77~^;2 2^2 +7~^ 2' (C^) 

where np is the Fermi occupation factor. It is straightforward to check that the exact Green function and the "raw" 
approximation given by the above value of F 

G« (*^) = 1^ (C4) 

coincide up to order to included. But this expansion of F, although undoubtedly correct, leads to the general causality 
problem described in Sect. Ill A. On the other hand, the reconstructed continued fraction is 

gS'^M^ -2 . (C5) 

+ (1 - (3unp{u))to - 

Fig.O shows a comparison between the spectral weight of the exact solution, the raw approximation, and the final 
continued fraction. One can see that Gj (iw) yields an extremely good spectral function (which remains true at any 
temperature and for a very wide range of values of to). 

Nonetheless, our enthusiasm has to be tempered by the following remarks. First, the exact solution being a two-pole 
rational function, it is not surprising that a finite continued fraction is able to mimic it satisfactorily. Secondly, going 
to low temperature does no harm in this case, but the two-site problem, described in Sect. D 2, shows that this is not 
true in general. 
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APPENDIX D: HIGHER LOOPS AND THE TWO-SITE PROBLEM 



1. Order t"^ and order 



Fig. |T8| shows the diagrams contributing to fourth and fifth order. The diagrams of Fig. |l^ lead to the foUowing 
expression for F: 



+ t 



6dfi 



3P'^d{-l + 2d) e'^iiuj) 



U(l + eO (M^-l)^(M^-l)^ 

3/3d(-l + e^) iuj 



( (l + e'^f{iu; - ifiiu; + if (l + e/f{iiu - lf{tiu + if 
X (-2 - 13e^ + 18rfe^ - 2e^'^ + 2iiu;f + e^{iujf + Me^iiujf + 2e^^{iujf) 

+ 2 1 1 

2(1 + e'^) (iuj - 3) {iu - lf{iu + if {iu + 3) 

X (247 - 278d + OUe'^ - 1396de'^ + 2476^'^ - 278^6^'^) 

+ (15 - 378d + 450e'^ - 1596de>^ + I5e^>^ - 378^6^^) {icuf 
+ (-3 + 42d - 90p/ + 252de^ - 3e'^'^ + A2de^'^){iujf 

+ (-3 + 6d + 6e'3 - I2de^ - Ze^<^ + Gde^^) {iujf 



(Dl) 



where we have set u = \. The fifth or|der and the corresponding continued fraction are too lengthy to be presented 
here, but are available on the internet.E3 
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FIG. 19. Above: Spectral function A{k = 0,tj) of the two-site Hubbard model with hopping term t = 0.8u at temperature 
T = t = 0.8u. Comparison between the exact result (bold) and the continued fraction at third (solid), fourth (dashed) and fifth 
(dot-dashed) order. In this region of parameter space, the approximation improves with increasing order. Below: The same, 
but with t = 0.5u at temperature T = O.lt — 0.05u. In this region of parameter space, the third-order approximation is the 
best one. 



2. The two-site problem 



The exactly solvable two-site problem provides a good testing ground for the higher-order results just obtained. In 
this subsection the parameter t is replaced by t/2 in order to account for periodic boundary conditions on a lattice 
with only two sites (say site 1 and site 2). The analytic expressions Qii{iu!) and 5i2(iw) of the Green function are 
available. Replacing the momentum integrations by sums over the two values A: = and k = n, the expansion of Qn 
up to order has been computed in the usual way (similar to what was done in the previous subsection), and checked 
to be the same as the one obtained from the exact solution. Furthermore, it leads to a continued fraction having only 
four floors instead of eight, as it should since the exact solution has four poles. It is then possible to investigate the 
validity of the various approximations throughout the (T, t) plane. Empirically, we were led to following conclusion: 
for T > t the approximation manifestly improves with the order, and for T <t the approximation deteriorates much 
faster with increasing order. Two examples illustrating this behavior are shown in Fig. 



1 


L. 


2 


L. 


3 


J. 


4 


E. 


5 


S. 


6 


J. 


7 


J. 


8 


C. 


9 


R. 


10 


H. 


11 


H. 


12 


H. 


13 


A. 



D. Landau, Sov. Phys. JETP 3, 920 (1957). 
D. Landau, Sov. Phys. JETP 5, 101 (1957). 
Hubbard, Proc. R. Soc. (London) Ser. A 276, 238 (1963). 
H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968). 
Tomonaga, Prog. Theor. Phys. 5, 544 (1950). 
M. Luttinger, J. Math. Phys. 4, 1154 (1963). 
Solyom, Adv. Phys. 28, 201 (1979). 

Bourbonnais and L. G. Caron, Int. J. Mod. Phys. B 5, 1033 (1991). 
Shankar, Rev. Mod. Phys. 66, 129 (1994). 
Frahm and V. E. Korepin, Phys. Rev. B 42, 10553 (1990). 
Frahm and V. E. Korepin, Phys. Rev. B 43, 5653 (1991). 

J. Schulz, in Les Houches, Session LXI (1994), edited by E. Akkermans et al. (Elsevier, Amsterdam, 1995), p. 533. 
V. Chubukov, D. Pines, and B. P. Stojkovic, J. Phys. Cond. Matt. 8, 10017 (1996), and Refs. lye 1992 and Ong 1990 
included therein. 

F. Gebhard, The Mott Metal-Insulator Transition (Springer, Berlin, 1997). 
J. Favand et al, Phys. Rev. B 55, R4859 (1997). 

S. E. Barnes, J. Phys. F 6, 1375 (1976). 
P. Coleman, Phys. Rev. B 28, 5255 (1983). 

G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57, 1362 (1986). 
Y. M. Vilk and A.-M. S. Tremblay, J. Phys. I France 7, 1309 (1997). 
N. Bulut, D. J. Scalapino, and S. R. White, Phys. Rev. B 50, 7215 (1994). 
E. Dagotto, A. Nazarenko, and M. Boninsegni, Phys. Rev. Lett. 73, 728 (1994). 
K. Kobayashi et al., Phys. Rev. Lett. 82, 803 (1999). 
A. Moreo et al, Phys. Rev. B 41, 2313 (1990). 
R. Preuss et at, Phys. Rev. Lett. 73, 732 (1994). 

R. Preuss, W. Hanke, and W. von der Linden, Phys. Rev. Lett. 75, 1344 (1995). 
W. Metzner and D. VoUhardt, Phys. Rev. Lett. 62, 324 (1989). 
T. Pruschke, D. L. Cox, and M. Jarell, Phys. Rev. B 47, 3553 (1993). 
A. Georges et al. Rev. Mod. Phys. 68, 13 (1996). 

D. E. Logan and P. Nozieres, Phil. Trans. R. Soc. Lond. A 356, 249 (1998). 
R. M. Noack and F. Gebhard, Phys. Rev. Lett. 82, 1915 (1999). 
K. Kubo, Prog. Theor. Phys. 64, 758 (1980). 
^2 K. K. Pan and Y. L. Wang, Phys. Rev. B 43, 3706 (1991). 
M. Bartkowiak and K. A. Chao, Phys. Rev. B 46, 9228 (1992). 



27 



C. Bourbonnais, Ph.D. thesis, Universite de Sherbrooke, 1985. 
S. K. Sarker, J. Phys. C 21, L667 (1988). 

S. Pairault, D. Senechal, and A.-M. S. Tremblay, Phys. Rev. Lett. 80, 5389 (1998). 
J. Voit, Rep. Prog. Phys. 57, 977 (1994). 

J. Voit, in Proceedings of The Ninth International Conference on Recent Progress in Many-Body Theories, Sydney, edited 
by D. Neilson (World Scientific, Singapore, 1998). 

C. Kim et at, Phys. Rev. Lett. 77, 4054 (1996). 
K. Kobayashi et at, Phys. Rev. Lett. 82, 803 (1999). 
B. O. Wells et at, Phys. Rev. Lett. 74, 964 (1995). 

P. W. Leung, B. O. Wells, and R. J. Gooding, Phys. Rev. B 56, 6320 (1997). 
W. Metzner, Phys. Rev. B 43, 8549 (1991). 

D. Boies, C. Bourbonnais, and A.-M. S. Tremblay, Phys. Rev. Lett. 74, 968 (1995). 
The 1/(7?!)'^ prefactor corrects Eq. (5) of Ref. ^ where l/(n!)^ should have appeared. 
J. Gilewicz, Approximants de Pade (Springer- Verlag, Berlin, 1978), Vol. 667. 
A. Georges and W. Krauth, Phys. Rev. B 48, 7167 (1993). 

J. W. Negele and H. Orland, Quantum many-particle systems (Addison- Wesley, New York, 1988). 
D. S. D. Poulin, S. Pairault and A.-M. S. Tremblay (unpubhshed). 
M. Jarrell and J. Gubernatis, Physics Reports 269, 133 (1996). 
A. Georges and W. Krauth, Phys. Rev. Lett. 69, 1240 (1992). 
S. Allen, private communication (unpublished). 
A.-M. Dare and G. Albinet, private communication (unpublished). 

V. S. Viswanath and G. Miiller, The Recursion Method, Lecture Notes m Physics (Springer- Verlag, Berlin, New York, 1994). 
®^ A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Dover, 
New York, 1963). 

S. Pairault, Mathematical" notebook at littp : //www. physique .usherb . ca/~dsenech/pairault .nb. 



28 



